package scu.maqiang.elastodynamics;

public class PEDynamicsSOTS1DL2 {

	public static void main(String[] args) {
//		long begin = System.nanoTime();
//		int N = 8;
//		double ep = 1.0 / N;
//
//		ScalarFunc coefG1 = (xyz, label, param) -> label == 10? 200: 4;
//		ScalarFunc coefEp1 = (xyz, label, param) -> label == 10? 1.0: 0.02;
//		ScalarFunc coefE1 = (xyz, label, param) -> label == 10? 5.0: 0.3;
//		ScalarFunc coefrho = (xyz, label, param) -> label == 10? 1: 0.02;
//		ScalarFunc coeff1 = (xyz, label, param) -> label == 10? 1: 1;
//		ScalarFunc coefrf = (xyz, label, param) -> label == 10? 1: 1;
//		//���㵥������
//		Mesh1L2Generator mg = new Mesh1L2Generator();
//		Mesh1DL2 cellMesh = mg.line(0, 1.0, 256);
//		ScalarFunc cellRegionFunc1D = (x, llable, param) -> {
//			return Math.abs(x[0] - 0.5) < 0.25?1.0:-1.0;
//		};
//		cellMesh.setDomainLabel(cellRegionFunc1D, null, 10);
//
//		FES1L2PE fsCell = new FES1L2PE(cellMesh);
//		double[][] FOCS = new double[2][fsCell.GetNdof()];
//		double[][] SOCS = new double[3][fsCell.GetNdof()];
//		double[][] D = new double[2][3];
//		SOTS1DL2PE ss = new SOTS1DL2PE(fsCell);
//		ss.PECellSolution(new ScalarFunc[] { coefG1, coefEp1, coefE1, coefrho, coeff1, coefrf }, null, FOCS, D, SOCS);
//		System.out.println("OK");
//		double time = 0.0, dt = 0.005, endTime = 10.0;
//		double alpha = 0.25, delta = 0.5;
//		double c0 = 1.0 / (alpha * dt * dt);
//		double c1 = delta / (alpha * dt);
//		double c2 = 1.00 / (alpha * dt);
//		double c3 = 0.50 / alpha - 1.0;
//		double c4 = delta / alpha - 1.0;
//		double c5 = 0.5 * dt * (delta / alpha - 2.0);
//		double c6 = dt * (1 - delta);
//		double c7 = delta * dt;
//
//		Mesh1DL2 homoMesh = mg.line(0, 1, 64);
//		FES1L2PE fsHomo = new FES1L2PE(homoMesh);
//		SRMatrix KHomo = new SRMatrix(fsHomo.GetNdof());
//		SRMatrix M0Homo = new SRMatrix(fsHomo.GetNdof());
//		SRMatrix MHomo = new SRMatrix(fsHomo.GetNdof());
//		double[] RHSHomo = new double[fsHomo.GetNdof()];
//		double[] uHomo = new double[fsHomo.GetNdof()];
//		double[] vHomo = uHomo.clone();
//		double[] aHomo = uHomo.clone();
//		double[] u0Homo = uHomo.clone();
//		double[] v0Homo = uHomo.clone();
//		double[] a0Homo = uHomo.clone();
//
//		fsHomo.assemblePEStiff(new double[] { D[0][0], D[0][1], D[0][2] }, BVPType.COMMON, KHomo);
//		fsHomo.assemblePEMass(new double[] {D[1][0]}, BVPType.CONSIST_MATRIX, M0Homo);
//		fsHomo.assemblePEMass(new double[] {D[1][0]}, BVPType.CONSIST_MATRIX, MHomo);
//
//		SRMatrix LHomo = KHomo.add(c0, MHomo);
//		double[] uvaHomo = uHomo.clone();
//
//		ScalarFunc initFunc =  (x, label, param) -> Math.sin(0.5 * Math.PI * x[0]);
//		double[] uu0Homo = fsHomo.valueFromFunc(initFunc, null);
//
//		for(int i = 0; i < uu0Homo.length; i++) {
//			u0Homo[2 * i] = uu0Homo[i];
//		}
//		double[] Ku0Homo = KHomo.mul(u0Homo);
//
//		DirectSSolver mHomoSolver = new DirectSSolver(M0Homo, fsHomo.getnBand());
//		//M0.sort();
//		//System.out.println(M0);
//		//System.exit(0);
//		for(int i = 0; i < fsHomo.GetNdof(); i++) {
//			RHSHomo[i] = RHSHomo[i] - Ku0Homo[i];
//		}
//
//		mHomoSolver.PLUSolve(RHSHomo, a0Homo, 1);
//		Arrays.setAll(aHomo, i -> a0Homo[i]);
//
//		for(int i = 0; i < fsHomo.GetNdof(); i++) {
//			uvaHomo[i] = c0 * u0Homo[i] + c2 * v0Homo[i] + c3 * a0Homo[i];
//		}
//		MHomo.mul(uvaHomo, RHSHomo);
////		String outputFile1 = "PEFreeVibrationSOTS1Du.dat";
////		String outputFile2 = "PEFreeVibrationSOTS1Dv.dat";
////		time = 0.0;
////
////		double[][] uvtimeHomo = new double[2][homoMesh.getNv()];
////		fsHomo.extractComponent(u0Homo, uvtimeHomo);
////		homoMesh.toTecplot(outputFile1, time, uvtimeHomo[0]);
////		homoMesh.toTecplot(outputFile2, time, uvtimeHomo[1]);
//		fsHomo.applyBC_MBN(LHomo, Direct.X, 1);
//		fsHomo.applyBC_MBN(LHomo, Direct.Y, 2);
//
//		DirectSSolver homoSolver = new DirectSSolver(LHomo, fsHomo.getnBand());
//
//		Mesh1DL2 fineMesh = mg.line(0, 1, 64);
//		ScalarFunc regionFineFunc1D = (x, llable, param) -> {
//			double xCell = x[0] * N - (int) (x[0] * N);
//			return Math.abs(xCell - 0.5) < 0.25 ? 1.0 : -1.0;
//		};
//		fineMesh.setDomainLabel(regionFineFunc1D, null, 10);
//
//		FES1L2PE fsFine = new FES1L2PE(fineMesh);
//		SRMatrix KFine = new SRMatrix(fsFine.GetNdof());
//		SRMatrix M0Fine = new SRMatrix(fsFine.GetNdof());
//		SRMatrix MFine = new SRMatrix(fsFine.GetNdof());
//		double[] RHSFine = new double[fsFine.GetNdof()];
//		double[] uFine = new double[fsFine.GetNdof()];
//		double[] vFine = uFine.clone();
//		double[] aFine = uFine.clone();
//		double[] u0Fine = uFine.clone();
//		double[] v0Fine = uFine.clone();
//		double[] a0Fine = uFine.clone();
//
//		fsFine.assemblePEStiff(new ScalarFunc[] { coefG1, coefEp1, coefE1 }, null, BVPType.COMMON, KFine);
//		fsFine.assemblePEMass(new ScalarFunc[] {coefrho}, null, BVPType.LUMP_MATRIX, M0Fine);
//		fsFine.assemblePEMass(new ScalarFunc[] {coefrho}, null, BVPType.CONSIST_MATRIX, MFine);
//
//		SRMatrix LFine = KFine.add(c0, MFine);
//		double[] uvaFine = uFine.clone();
//
//		double[] uu0Fine = fsFine.valueFromFunc(initFunc, null);
//
//		for(int i = 0; i < uu0Fine.length; i++) {
//			u0Fine[2 * i] = uu0Fine[i];
//		}
//		double[] Ku0Fine = KFine.mul(u0Fine);
//
//		DirectSSolver mFineSolver = new DirectSSolver(M0Fine, fsFine.getnBand());
//
//		for(int i = 0; i < fsFine.GetNdof(); i++) {
//			RHSFine[i] = RHSFine[i] - Ku0Fine[i];
//		}
//
//
//		System.out.println(Arrays.toString(RHSHomo));
//		System.out.println(Arrays.toString(RHSFine));
//		System.exit(0);
//		mFineSolver.PLUSolve(RHSFine, a0Fine, 1);
//		Arrays.setAll(aFine, i -> a0Fine[i]);
//
//		for(int i = 0; i < fsFine.GetNdof(); i++) {
//			uvaFine[i] = c0 * u0Fine[i] + c2 * v0Fine[i] + c3 * a0Fine[i];
//		}
//		MFine.mul(uvaFine, RHSFine);
//		time = 0.0;
//		fsFine.applyBC_MBN(LFine, Direct.X, 1);
//		fsFine.applyBC_MBN(LFine, Direct.Y, 2);
//
//		DirectSSolver fineSolver = new DirectSSolver(LFine, fsFine.getnBand());
//
//		double[] fineU0 = new double[fsFine.GetNdof()];
//		double[] fineU1 = new double[fsFine.GetNdof()];
//		double[] fineU2 = new double[fsFine.GetNdof()];
//
//		double[][] fineU_uphi = new double[2][fineMesh.getNv()];
//		double[][] fineU0_uphi = new double[2][fineMesh.getNv()];
//		double[][] fineU1_uphi = new double[2][fineMesh.getNv()];
//		double[][] fineU2_uphi = new double[2][fineMesh.getNv()];
//		double[][] u012 = new double[4][];
//		double[][] phi012 = new double[4][];
//
//		String outputFile1 = "PEDynamicsSOTS1Du.dat";
//		String outputFile2 = "PEDynamicsSOTS1Dphi.dat";
//
//		fsFine.extractComponent(uFine, fineU_uphi);
//		u012[0] = fineU_uphi[0];
//		u012[1] = fineU_uphi[0];
//		u012[2] = fineU_uphi[0];
//		u012[3] = fineU_uphi[0];
//
//		phi012[0] = fineU_uphi[1];
//		phi012[1] = fineU_uphi[1];
//		phi012[2] = fineU_uphi[1];
//		phi012[3] = fineU_uphi[1];
//		fineMesh.toTecplot(outputFile1, time, u012);
//		fineMesh.toTecplot(outputFile2, time, phi012);
//		double[] myParam = new double[1];
//		while(time <= endTime) {
//			time += dt;
//			myParam[0] = time;
//			System.out.println("time: " + time);
//			//fs.applyBC_MBN(RHS, Direct.All, 0.0, 1, 2);
//			fsHomo.applyBC_MBN(RHSHomo, Direct.X, 0.0, 1);
//			fsHomo.applyBC_MBN(RHSHomo, Direct.Y, 0.0, 2);
//			fsFine.applyBC_MBN(RHSFine, Direct.X, 0.0, 1);
//			fsFine.applyBC_MBN(RHSFine, Direct.Y, 0.0, 2);
//			//System.out.println(Arrays.toString(RHS));
//
//			homoSolver.PLUSolve(RHSHomo, uHomo, 1);
//			fineSolver.PLUSolve(RHSFine, uFine, 1);
//			for(int i = 0; i < fsHomo.GetNdof(); i++) {
//			    aHomo[i] = c0 * (uHomo[i] - u0Homo[i]) - c2 * v0Homo[i] - c3 * a0Homo[i];
//			    vHomo[i] = v0Homo[i] + c6 * a0Homo[i] + c7 * aHomo[i];
//			    uvaHomo[i] = c0 * uHomo[i] + c2 * vHomo[i] + c3 * aHomo[i];
//			}
//
//			for(int i = 0; i < fsFine.GetNdof(); i++) {
//			    aFine[i] = c0 * (uFine[i] - u0Fine[i]) - c2 * v0Fine[i] - c3 * a0Fine[i];
//			    vFine[i] = v0Fine[i] + c6 * a0Fine[i] + c7 * aFine[i];
//			    uvaFine[i] = c0 * uFine[i] + c2 * vFine[i] + c3 * aFine[i];
//			}
//
//			MHomo.mul(uvaHomo, RHSHomo);
//			MFine.mul(uvaFine, RHSFine);
//			for(int i = 0; i < fsHomo.GetNdof(); i++) {
//				u0Homo[i] = uHomo[i];
//			    v0Homo[i] = vHomo[i];
//			    a0Homo[i] = aHomo[i];
//			}
//
//			for(int i = 0; i < fsFine.GetNdof(); i++) {
//				u0Fine[i] = uFine[i];
//			    v0Fine[i] = vFine[i];
//			    a0Fine[i] = aFine[i];
//			}
//
//
//			ss.assemblePESOTSSolution(fsFine, fineU0, fineU1, fineU2, fsHomo, uHomo, fsCell, FOCS, SOCS, N, ep);
//			fsFine.extractComponent(uFine, fineU_uphi);
//			fsFine.extractComponent(fineU0, fineU0_uphi);
//			fsFine.extractComponent(fineU1, fineU1_uphi);
//			fsFine.extractComponent(fineU2, fineU2_uphi);
//			u012[0] = fineU_uphi[0];
//			u012[1] = fineU0_uphi[0];
//			u012[2] = fineU1_uphi[0];
//			u012[3] = fineU2_uphi[0];
//
//			phi012[0] = fineU_uphi[1];
//			phi012[1] = fineU0_uphi[1];
//			phi012[2] = fineU1_uphi[1];
//			phi012[3] = fineU2_uphi[1];
//
//			fineMesh.toTecplot(outputFile1, time, u012);
//			fineMesh.toTecplot(outputFile2, time, phi012);
//		}
//
//		long end = System.nanoTime();
//		System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
	}

}
